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Abstract 

We present results for the gluon and ghost propagators in SU(N) Yang-Mills theory 
on a four-torus at zero and nonzero temperatures from a truncated set of Dyson- 
Schwinger equations. When compared to continuum solutions at zero temperature 
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hand, are minute for T < 250 MeV. 
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1 Introduction 



The infrared behaviour of QCD Green's functions is intimately related to 
the issue of confinement, see e.g. ref. [1]. The confinement phenomenon is 
intrinsically non-perturbative, and thus a quite limited number of methods is 
suited to investigate this and related issues. Furthermore, the infrared singular 
nature of some of the QCD Green's functions provides a further severe obstacle 
in the study of these quantities. In ref. [2] a numerical method is described 
where the analytical treatment of the gluon and ghost propagator Dyson- 
Schwinger equations (DSEs) for momenta smaller than a certain matching 
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scale allowed for the numerical solution of these coupled non-linear integral 
equations. At this point the question arises whether a purely numerical method 
to solve these equations alleviates their treatment. To this end one notes that 
on a compact manifold the finite volume of underlying space-time acts as a 
natural infrared regulator. Using a four-torus this will, on the one hand, be 
exploited here for a direct solution of the coupled gluon and ghost DSEs, 
and, on the other hand, be utilized to estimate effects of the finite volume. 
An extension to asymmetric tori is equivalent to the choice of non-vanishing 
temperatures, and thus it will also be considered in this paper. 

As noted above confinement manifests itself in the long range behaviour of 
Yang-Mills theory. According to the Kugo-Ojima confinement criterion [3,4] 
as well as the Gribov-Zwanziger horizon scenario, see e.g. ref. [5], it is expected 
that the gauge fixing degrees of freedom and thus the Faddeev-Popov ghosts 
provide the long-range correlations. In Landau gauge this has been confirmed 
in the framework of DSEs as well as exact renormalisation group equations 
in a series of papers [5,6,7,8,9,10,11,12]: the ghost propagator is indeed more 
singular in the infrared than a simple pole. The gluon propagator, on the 
other hand, vanishes in the infrared. Lattice Monte-Carlo calculations agree 
with these findings, although it is currently under debate whether the lattice 
gluon propagator is finite [13] or vanishes in the infrared [14]. Here a study of 
DSEs at finite volumes can be helpful when judging the finite-volume effects 
in lattice calculations. As due to the periodic boundary conditions used in 
lattice calculations the underlying space-time manifolds are effectively four- 
tori, we will employ these manifolds in our study. In this paper we will present 
general techniques needed to solve the related DSEs, and we will comment 
on the differences in the method and for the solutions when compared to the 
treatment based on flat Euclidean space-time [2,7,15,16]. 

The paper is organised as follows: In the next section we summarise general 
properties of DSEs on a compact manifold. We discuss the boundary condi- 
tions of the fields, present the equations for finite temperature and discuss the 
renormalisation. Numerical results are presented in section 3: Finite volume 
effects are discussed in subsection 3.1, and results for non- vanishing temper- 
atures are given in section 3.2. We conclude with a summary and an outlook 
on further applications in section 4. The employed integral kernels are listed 
in appendix A and details of the numerical method are presented in appendix 
B. 



2 Dyson-Schwinger equations on a compact manifold 

The DSEs are equations for the Green's functions, i.e. for the vacuum ex- 
pectation values of time-ordered products of quantum fields. Therefore it is 
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quite illustrative to consider first the general properties of quantum fields on 
compact space-time manifolds. 



2. 1 Quantum fields on a compact Euclidean space-time manifold 



Quantum fields on a compact Euclidean space-time manifold are realised by 
imposing ( ant i-) periodic boundary conditions: 

+ L M ) = e i2 ^$(^) . (1) 

Here <J> is shorthand for all fields that appear in the Lagrangian of the theory, 
and L M denotes the length of the torus in the /i-th direction. The overall 
phase is i] = for bosons, i.e. periodic boundary conditions, and r] — \ for 
fermions, i. e. antiperiodic boundary conditions. Furthermore, in <i-dimensional 
coordinate space <3>(x) has the following Fourier representation 

= T^—r E «V.~n d e^ , (2) 

Li\ - • • Lid ni...n d 

= — (% + ?7), (3) 

in terms of its Fourier components $ ni ... nd - The oj^ are discretized momenta, 
and due to an obvious analogy we will call these Matsubara momenta or 
Matsubara frequencies. The special case we want to study here is d = 4, i.e. 
three space dimensions (// = 1, 2, 3) and one 'time' (// = 4) dimension. 

In the framework of the imaginary time formalism quantum field theory at 
non-vanishing temperatures is related to the corresponding Euclidean path- 
integral with compactified time-direction such that At = (3 = 1/T. Therefore 
torus-compactified Euclidean space-time with a shorter length for the time 
direction introduces an effective temperature which can be estimated to be 

T= l/L 4 withL 4 <L . (4) 

Depending on the application we will either set all L M equal, or we choose 
L t = L 2 = L 3 = L and L 4 = (3 = 1/T. 

For flat Euclidean space-time the DSEs in momentum space are integral equa- 
tions. Due to the representation of the fields in Fourier space of compact 
space-time we deal in this case with sums over Matsubara momenta, i. e. when 
compactifying space-time to a torus the DSE integrals are replaced by sums 
over the n„ 
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Fig. 1. Sketch of the momentum grid dual to the four-torus. The summation over 
complete hyperspheres is indicated by fully drawn circles. The hyperspheres depicted 
by dashed lines are not complete if one uses Cartesian cutoffs instead of an 0(3) 
invariant one. Note that an 0(4)-invariant cutoff can only be chosen in the zero 
temperature case. 

For the numerical treatment of the equations it is convenient to rearrange this 
summation such that they represent a spherical coordinate system [17], see 
fig. 1 for an illustration. Employing the 0(3) symmetry of the space directions 
we write 



with n = 714. Furthermore, j denotes spheres with = const, i — 1 . . . 3, 
and m numbers the grid points on a given sphere. This spherical summation 
of Cartesian grid points is very useful device when introducing 0(3) invariant 
cut-offs. 

On a torus the gluon as a boson clearly obeys periodic boundary conditions. 
The ghost field occurs when representing the Faddeev-Popov determinant as a 
Grassmannian functional integral. Although the ghost fields c and c are thus 
anticommuting fields they also have to obey periodic boundary conditions. 
Including fermionic fields \1/ in the fundamental representation this can be 
seen easily from the ghost fields' properties under BRST-transformations: 



Here the t a are the generators of SU(N), and s is the Slavnov operator. Since 
s is continuous and thus unique, it is periodic. To preserve the antiperiodicity 
(due to the fermionic nature of the matter field) on both sides of the equation 
periodic boundary conditions for the ghost field c are required. 





stf = -igt a c a ^ . 



(7) 



2.2 The Dyson-Schwinger Equations 



The fully renormalized zero temperature ghost and gluon propagators in Lan- 
dau gauge can always be written in the form 



iJ„(t)=(^-^)^. 0) 

where G(/c 2 ) and Z(k 2 ) denote the ghost and gluon dressing functions, respec- 
tively. A combination of these functions can be used to define the nonpertur- 
bative running coupling [6] 

a(k 2 ) = a(/i 2 ) G 2 (k 2 , /i 2 ) Z(k 2 , /i 2 ). (10) 

which is a renormalization group invariant, i.e. independent of the renormal- 
ization scale //. 

For non-vanishing temperatures, in the rest frame of the heat-bath, one has 
two different tensor structures in the gluon propagator, one transverse and 
one longitudinal to the heat bath [18,19]: 



D tM = fcif Z m(k : \k\) + P^M Z (k , \k\)) With (11) 

P, L M = P,Ak)-P^k), P»» = ^-^, (12) 
i,3 = 1,2,3; fi,u = 1,2,3,4, 

and k 2 = k 2 , + k 2 . The ghost propagator is a Lorentz scalar and does not 
acquire further structures but its dressing function will depend then seperately 
on frequency and three-momentum: 

Sab 

B*(k) = -p-G(*b, \k\). (13) 
The following linear combination of gluon dressing functions, 

Z T (k):=±Z (k) + lz m (k), (14) 
will prove to be closest to the zero temperature gluon dressing function. 
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Fig. 2. Diagrammatic representation of the propagator DSEs in the truncation 
scheme used in this work. Wiggly lines denote heat-bath transverse gluon prop- 
agators, double-dashed lines are heat-bath longitudinal propagators and dashed 
lines represent ghost propagators. Blobs indicate dressed propagators and vertex 
functions, respectively. 

The diagrammatic structure of the Dyson-Schwinger equations for the ghost 
and gluon propagators is shown in fig. 2. It is a coupled, nonlinear system 
of equations, which contains dressed propagators as well as dressed vertex 
functions on the right hand side. In order to obtain a closed system of equations 
it is necessary to specify suitable approximations for these vertices. Such a 
truncation scheme has been developed and successfully applied in ref. [17]. As 
the scheme is discussed in great detail in ref. [15], we only mention briefly 
its main ideas here. The key observation is ghost dominance in the infrared 
[12], i.e. the ghost loop in the gluon DSE dominates all contributions from 
the gluon diagrams for small momenta. It is furthermore known, that the bare 
ghost-gluon vertex, 

rf- 9 \p,q)=tq» (15) 
is an excellent approximation in both the infrared and ultraviolet momentum 
regime [20,21]. Here q^ is the momentum of the outgoing ghost-leg and the 
colour-indices have been suppressed. The gluon two-loop diagrams, involving 
the full four-gluon vertex, can be neglected in these regions and it turned 
out to be a suitable approximations to neglect contributions from the four- 
gluon vertex altogether. For the three-gluon vertex a minimal dressing has 
been chosen such that the correct perturbative anomalous dimensions of the 
propagators in the ultraviolet are recovered in the solutions. The vertex is 
given by 

r %\( k ,P,q) = T3 J 3 fx(k,p,q)H 3g ( y k,p,q) (16) 
where T 3 ^(k,p, q) is the bare vertex. The dressing function H Sg is given by 

U 3g {k,p,q) - — z{q){1+35) z(p) (i +35 ) > ( 17 ) 
with the momenta q and p running inside the loop. Here the anomalous dimen- 
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sion of the ghost propagator, 8 = —9/44, together with the vertex renormaliza- 
tion constant Z\ ensure the correct ultraviolet running of the vertex with mo- 
menta and renormalization scale. Furthermore this choice leads to cutoff inde- 
pendent ghost and gluon dressing functions in the continuum. This truncation 
scheme has been generalised to non- vanishing temperatures in ref. [22,23,24]. 

On an asymmetric torus this leads to the following set of DSEs for the ghost 
and gluon dressing function G(k), Z m (k) and Z (k): 

W)=* 3+ g2NcZi ^ £ m w^h? {ATZm{k ~ q) + AhUk ~ q)) ' 

(18) 

1 7 1 n 2 N 7 T V" p i 1 n 2 M 7 T V" H 3g {q,p,k) 

— -ZS--S N cZl - £ - 1 ^ r R + -g N cZl - £ 

(M T Z m (q)Z m (p) + M 1 Z (q)Z m (p) + M 2 Z (p)Z m (q) + M L Z (q)Z (p)) , 

(19) 

(N T Z m (q)Z m {p) + N 1 Z (q)Z m {p) + N 2 Z (p)Z m (q) + N L Z (q)Z (p)) . 

(20) 

The expressions for the kernel functions At, Al, R, Mt, M\, M 2 Ml, P and 
Nt, Ni, N 2 , Nl are given in appendix A, see eqs. (A.1)-(A.12). The corre- 
sponding 0(4)-symmetric version of these equations in the zero temperature 
limit can be found in ref. [17]. 



2.3 Renormalisation 



We employ a momentum subtraction (MOM) scheme to renormalise the DSEs 
on the torus. For zero temperature, i.e. on a symmetric torus, this scheme has 
been described in [17]. The DSEs for the (renormalised) dressing functions 
D G {G, Z m , Z } at the four momentum k 2 can be written symbolically as 

^p^ + n^), (2i) 

where Z D stands for the respective renormalisation constants Z 3 , Z 3m , Z 30 
and Ur){k 2 ) denote the loop contributions. The momenta kf are lying on the 
hyperspheres shown in fig. 1. To eliminate the momentum independent renor- 
malisation factors Z D we subtract eqs. (21) from their respective equations 
evaluated at a fixed subtraction scale s. This scale is conveniently chosen equal 
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to the squared renormalisation point, s = fi 2 . We obtain: 



1 1 



+ U D (k 2 )-U D ^ 2 ) . 



(22) 



D(k 2 ) D(fi 2 ) 



One then has to specify a renormalisation condition, i.e. a value for D(fj, 2 ), 
to complete the (re-)normalisation procedure. For convenience this condi- 
tion has been adapted from the continuum solution, given in ref. [10]. Using 
a(/j, 2 ) = 0.968 the renormalisation point is given by fi 2 = 1.610 GeV 2 and 
the normalisation condition is Z(fi 2 ) = 0.847. The corresponding value of the 
ghost dressing function is then G(/i 2 ) = l/^jz(fx 2 ) = 1.087, c.f. eq. (10). 

On the asymmetric torus we exploit the fact that no new divergences can 
arise at finite temperatures in a renormalisable quantum field theory [18,19]. 
Only the finite part of the renormalisation, i.e. the condition on D(s), can be 
affected. If the subtraction scale s is much larger than maxjAqcD, 1/T}, where 
the dressing functions are not affected by temperature effects, it should suffice 
to determine only one renormalisation constant for all Matsubara frequencies 
of each respective propagator. These constants are determined e.g. from the 
zeroth order term in the Matsubara sum (n = 0) by the condition 



where D(0,s) can be conveniently set to the zero temperature value D(s). 
We checked that within a very small error we obtain the same value for the 
constant Z D also for all other Matsubara frequencies. The overall temperature 
dependence of the renormalisation constants is several percent for tempera- 
tures in the range from 50 MeV to 800 MeV. The renormalisation constants 
for the heat-bath transverse and longitudinal gluon dressing functions differ 
slightly, but at most a few percent for increasing temperatures. 

In regularisation schemes with an ultraviolet cutoff one also has to deal with 
quadratic divergencies, which cannot be accounted for by the MOM-procedure 
described above. In the continuum formulation of the DSEs quadratic diver- 
gencies are removed on the level of the DSEs by a careful procedure [10,22] 
which amounts to adding a suitable counterterm to the bare action of the the- 
ory. A similar method has been applied to the formulation on the symmetric 
and asymmetric torus and is described in detail in [17,23]. We therefore refrain 
from repeating the details here and just note that the integral kernels given 
in appendix A contain kernels where the dangerous terms have already been 
eliminated. 

The dressing functions are calculated on a momentum grid that is limited 
by the ultraviolet cut-off A. However, as can be seen directly in equations 
(18)-(20), the integral kernels and the dressing functions in the loops depend 
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p 2 [GeV 2 ] 

Fig. 3. The gluon propagator Z(p 2 )/p 2 in the continuum and on the torus at two 
different volumes. 



A[GeV] 


2.47 


3.09 


3.71 


4.33 


4.95 


Z)(0)[GeV- 2 ] 


4.90 


4.55 


4.23 


3.97 


3.75 



Table 1 " ' 

Finite size analysis: Dependence of -D(O) on the cutoff A at a fixed volume of V = 
4091 fm 4 . 



on the relative momentum between the momentum k from outside the loop 
and the loop momentum q. These relative momenta may have values between 
zero and two times the momentum cut-off. Thus one has to extrapolate the 
dressing functions beyond the cut-off. This is accomplished by employing the 
ultraviolet asymptotic forms (B.12) and (B.13), which are given in appendix 
B.2. These forms are known to be analytic solutions of the continuum DSEs 
at large momenta. For typical torus sizes one does not exactly reproduce the 
one-loop scaling at large momenta, since the momentum cutoff is still rather 
low of the order of a few GeV. However, we need the extrapolation only for 
a small number of momentum configurations and our procedure therefore is 
satisfactory. We also checked other prescriptions, like setting the dressing func- 
tions to a constant value or even to zero for momenta larger than the cut-off 
and found only very minor changes in the solutions of the DSEs. 
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y[fm 4 ] 


105964 


51102 


20931 


4091 


2461 


G(0) 


3.84 


3.79 


3.68 


3.36 


3.10 



Table 2 " " 

Finite volume analysis: Dependence of G(0) on the volume V at a fixed cutoff of 
A = 2.47 GeV. 

3 Results 

3.1 Finite volume effects at zero temperature 

Our first application of the numerical tools described in the previous sections 
is to determine whether one can recover the solutions of the continuum DSEs, 
determined in ref. [10], also on a torus. First results of such an investigation 
have already been reported in refs. [10,25]. Here we extend this investigation 
to much larger volumes seeking for signals of the continuum limit on the 
torus. To clarify the effects in the most pronounced way we display results 
for the gluon propagator, D(p 2 ) = Z(p 2 )/p 2 , in fig. 3 and the ghost dressing 
function, G(p 2 ), in fig. 4. The momentum scale on the torus is adapted to 
the one from the corresponding continuum calculation [17]. Our strategy was 
to fix the cutoff at a certain scale, A = 2.47 GeV for the solutions presented 
in figures 3 and 4, and then vary the numbers of points on the momentum 
lattice to generate different volumes of the torus. For p in units of GeV and 
denoting the smallest momentum accessible on the (symmetric) lattice, the 
corresponding volume is determined by V = (2tc/p • 0.1973 fm) 4 . The two 
solutions displayed in figures 3 and 4 have been calculated with V ~ 2500 fm 4 
and V w 10 5 fm 4 . Note that the smaller volume corresponds to a typical one 
achieved in recent lattice simulations [13], whereas the larger value exceeds the 
current possibilies of the lattice formulation by far. Interestingly, our results for 
these two volumes are very similiar to each other, whereas there is a qualitative 
gap to the continuum results. The numerical continuum solution for the gluon 
propagator is in accordance with the infrared analytical result, eq. (B.10), 
and clearly shows a vanishing gluon propagator at zero momentum with the 
anomalous dimension k ~ 0.595. The solutions on the torus, on the other 
hand, seem to go to a constant if extrapolated to zero momentum. 

In order to analyse the zero momentum limit more precisely we fitted the 
infrared part of our solution to the form 

Z ir ( P 2 )_aof P 2 V a \ 
p 2 p 2 \a 2 + p 2 J 

For momenta below p ~ 250 MeV we found very good agreement with our large 
volume solution for the values of the parameters a = 0.244, a 1 = 0.233 GeV 
and a2 = 0.500. Thus, to very good precision we find that the torus solution 
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p 2 [GeV 2 ] 

Fig. 4. The ghost dressing function G(p 2 ) in the continuum and on the torus at two 
different volumes. 

for the gluon propagator is constant in the infrared, and thus very different 
from its continuum limit. However, there is a caveat here. Finite size effects 
due to the ultraviolet momentum cutoff A uv could introduce contributions of 
the order 1/Auv or 1/A UV to D(0) on the torus. Provided these terms have 
large coefficients they could dominate even for very small momenta and hide 
a possible underlying power law with k > 0.5. In order to test for finite size 
effects we solved the DSEs on a torus with fixed volume V = 4091 fm 4 for five 
different cutoffs corresponding to momentum lattices of 16 4 , 20 4 , 24 4 , 28 4 and 
32 4 . We then fitted the expression (24) in the infrared with a 2 fixed to 0.5. 
Our results for D(0) are given in table 1. We find a considerable dependence 
of D(0) on the cutoff, which can be described by the form 

D(0) — do -\- ~~r~~ (25) 
Auv 

with the parameters do = 2.67GeV" 2 and d\ = 5.57G6V" 1 . If we use A uv 
instead of Auv the quality of the fit gets much worse, indicating that the 
dominant finite size effects are of order 1/Auv- Nevertheless, the nonvanishing, 
cutoff independent constant d indicates that the gluon propagator on a torus 
is indeed finite in the infrared. 

Our solutions for the ghost dressing function on a torus are shown in fig. 4. 
Similar to the gluon we note a clear qualitative difference between the torus 
and the continuum solutions. Although there is a distinct increase of G{p 2 ) 
in the midmomentum region, the function is much less steep in the infrared 
than the continuum solution. The volume dependence for a given cutoff is very 
small and it is not at all clear that G(0) — > oo for V — > oo. On the contrary, 
the infrared part of the ghost dressing function for the large volume can be 
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Fig. 5. The gluon propagator Z(p 2 )/p 2 in the continuum and on the torus compared 
to the lattice results of refs. [13,26,27]. The most recent lattice data from ref. [28] 
agree well with the data shown. 
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Fig. 6. The ghost dressing function G(p 2 ) in the continuum and on the torus com- 
pared to the lattice results of refs. [28,29]. 



fitted very well with the form 

G ir (p 2 ) 



(ei +p 2 ) e *' 



(26) 



where e = O.llGeV 262 , e x = 0.02GeV 2 and e 2 = 0.32. If the scale e x is 
not chosen as a fit parameter but set to zero the fit becomes unstable. This 
clearly indicates that the dressing function is finite when extrapolated to zero 
momentum. From this fit form we can determine G(0) for several volumes. 
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p 2 [GeV 2 ] 

Fig. 7. The running coupling in the continuum and on the torus compared to the 
lattice results of ref. [28]. 

Our results are given in table 2. From these values it is not clear whether the 
ghost diverges in the infinite volume limit. If so, then this limit is approached 
at an extremely slow rate. A possible polynomial representation of the values 
in table 2 is given by G(0) = go — g\/V 92 with g = 3.89, g\ = 0.305 and 
gi = 0.668, which indeed indicates a constant ghost in the infinite volume 
limit. However, other forms may fit the values equally well. 

In figs. 5 and 6 we compare our results for the gluon propagator and the 
ghost dressing function with corresponding ones from recent lattice calcula- 
tions [13,26,27,28,29] (see also [30] for first results on large asymmetric lat- 
tices). Interestingly enough, there is a striking qualitative agreement between 
the lattice data and our torus result for the gluon propagator. Both seem to go 
to a finite value in the infrared, in disagreement with the continuum solution. 
Thus, provided the analogy between the lattice formulation and the formula- 
tion of the DSEs on the torus holds, then the 'true', continuum limit gluon 
propagator may look quite different in the infrared than the one of contem- 
porary lattice calculations. As for the ghost, both results, the DSE solution 
in the continuum and the one on the torus look strikingly similar on a linear 
plot of the momentum range assessible on the lattice. In turn, this means that 
the lattice ghost propagator could equally well be divergent or finite at zero 
momentum. 

The running coupling (10), determined from a combination of the ghost and 
gluon dressing functions, is shown in fig. 7. Whereas in the continuum the 
coupling approaches the infrared fixed point (B.ll) for small momenta, we 
find an infrared vanishing coupling on the torus. Note that the difference 
between the small volume (12 4 ) and large volume (36 4 ) result on the torus is 
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Fig. 8. The zeroth Matsubara component of the gluon propagator, Zx/k^, in the 
infrared region at different temperatures (24 3 momentum grid). 




1.0 
k 3 [GeV] 



Fig. 9. The ghost dressing function G in the infrared region at different temperatures 
(24 3 momentum grid). 



very small. Thus one cannot extrapolate to the continuum limit by varying 
the torus volume. Compared to the results of recent lattice calculations [28] we 
again find qualitative agreement between the torus solutions and the lattice 
data 1 . 
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3.2 Ghost and gluon propagators at temperatures below T c 

A convenient way to assess the temperature dependence of the gluon dress- 
ing functions is to evaluate the linear combination Zx{k), eq. (14), and the 
difference of the dressing functions Z and Z m , 

AZ(k) = Z (k)-Z m (k). (27) 

The linear combination (14) connects directly to the gluon dressing function 
Z(p 2 ) at zero temperature whereas AZ measures temperature effects on the 
tensor structure of the gluon propagator. 

The obtained results have in common that the ghost and gluon propagators 
are only slightly temperature dependent for T < 240MeV. Hereby the most 
interesting parts are their zeroth Matsubara modes since all other modes turn 
out to be 0(4) invariant, i.e. they depend to a very high precision only on 
p 2 = p 2 + pi, and therefore show no explicit temperature effects. 

The results presented in figs. 8 and 9 have been calculated using a 24 3 grid 
and an ultraviolet cutoff of 4.7 GeV. This corresponds to a spatial volume of 
approximately 6 fm 3 . In figs. 8 and 9 we show the zeroth Matsubara modes of 
the ghost dressing function G and the quantity Zx/k 2 at several temperatures. 
As has been discussed in the previous section only the finite-volume value of 
the correspondig infrared exponents can be extracted. The propagator func- 
tion Zx/k 2 shows only small variations for a large range of temperatures up 
to at least 240 MeV. The ghost dressing functions also changes only by a few 
percent. The quantity AZ (27) which measures the change of the tensor struc- 
tur of the gluon propagator is displayed in fig. 10. Whereas AZ remains quite 
small for small temperatures, it increases drastically in the mid-momentum 
regime around A; 3 = 700 MeV when the temperature is raised above 200 MeV, 
thereby indicating a sizeable change of the gluon propagator in this region. 
These temperature effects are presumably most important just below the crit- 
ical temperature. Thus the question arises whether higher Matsubara modes 
are important for this quantity. As one sees from fig. 11 this is definitely not 
the case for T=140 MeV. At larger temperatures we find even smaller contri- 
butions from these modes. 

At T=140 MeV we have also studied the system for considerably larger vol- 
umes. As one can see from figs. 12 and 13 the changes when going to larger 
volumes are, even in the infrared, minute. Of course, this has to be expected 
from the results presented in previous subsection. Nevertheless, as discussed 
this does not mean that the infinite- volume limit of these quantities is reached. 

1 The quantitative difference between the two approaches at large momenta is due 
to missing two-loop contributions in the gluon-DSE. 
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Fig. 10. The difference AZ of the gluon dressing functions at different temperatures 
(24 3 momentum grid). 




Finally, we want to mention that with this numerical method it is possible to 
find qualitatively similar solutions for the gluon and ghost dressing functions 
also for quite high temperatures up to T=800 MeV [23]. At the moment we 
lack a definite interpretation of such solutions, the most likely explanation is 
that they belong to a super-heating of the low-temperature phase. 
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Fig. 12. The ghost dressing function G in the infrared region from different grid 
sizes at T=140 MeV). 




Fig. 13. The difference AZ of the gluon dressing functions from different grid sizes 
at T=140 MeV 

4 Summary 

We presented solutions for the coupled set of Dyson-Schwinger equations for 
the ghost and gluon propagators of Landau gauge Yang-Mills theory on a 
four-dimensional torus. We studied the influence of the finite volume on the 
infrared behaviour of the propagators on a symmetric torus and we employed 
asymmetric tori to introduce non-vanishing temperatures. 

Evaluating the ghost and gluon propagators on the torus for various vol- 
umes and momentum cutoffs provides a surprising result: Even for extremely 
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large volumes, exceeding contemporary lattice volumes by magnitudes, dis- 
tinct qualitative differences between torus and continuum solutions are found. 
On the torus, the gluon propagator as well as the ghost dressing function 
are finite in the infrared, whereas the same quantities in the continuum limit 
are infrared vanishing (gluon) or diverging (ghost). By varying the ultravio- 
let cutoff we made sure that these effects are not artefacts of the ultraviolet 
regularization. Thus the continuum limit seems to be a nontrivial transition 
from the compact to the noncompact manifold. 

We furthermore compared our results with contemporary lattice calculations. 
For the ghost dressing function all three solutions look similar in the mo- 
mentum range accessible by lattice calculations. For the gluon propagator, 
however, we found qualitative agreement of the torus results with the lat- 
tice data in the infrared, whereas the continuum result is different. Thus our 
results suggest that the continuum gluon propagator may be qualitatively 
different than the one found in contemporary lattice calculations. The same 
is true for the running coupling. The analytical results from the continuum 
DSE- approach give a fixed point in the infrared. Although the value of this 
fixed point is, within narrow limits, subject to the employed ansatz for the 
ghost-gluon vertex [9], its very existence can be derived from the structure of 
the ghost-DSE alone and is therefore independent of a truncation [31]. The 
fixed point is universal [12] and shown to be invariant in a class of gauges that 
interpolate between Landau and Coulomb gauge [32]. On a compact mani- 
fold, however, both, the DSE-approach and lattice Monte-Carlo simulations 
[28,33,34,35] find a vanishing coupling in the infrared in disagreement with 
the analytically derived continuum limit. 

A further problem studied in this work is the effect of non-vanishing temper- 
atures on the infrared anomalous dimensions of the ghost and gluon propaga- 
tors. Our results indicate that the anomalous dimensions are unaffected up to 
temperatures of the order of the phase transition. Moderate-size differences do 
occur in the tensor structures of the gluon propagator indicating a qualitative 
change in the gluon propagator as one approaches the phase transition from 
below. Nevertheless, it remains an open question which degrees of freedom are 
the most relevant for the deconfinement phase transition. 

We presented a numerical method to solve coupled sets of nonlinear integral 
equations on a compact manifold by iteration. The method is suitable for both 
symmetric and asymmetric space-time tori. It can be applied to systems with 
an arbitrary number of dimensions and an arbitrary number of equations. 
We are confident that the algorithm presented here also provides a sound 
basis for investigations of related issues. E.g. the choice of twisted boundary 
conditions on the torus [36] might hold the key to provide a tool for assessing 
the influence of topologically non-trivial field configurations on the infrared 
behaviour of QCD Green's functions. 
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A Integral kernels for the Dyson-Schwinger Equations 



The two kernels of the ghost equation (18) are given by: 



The kernel of the ghost loop and the four kernels of the gluon loop for the 
heat-bath transversal part of the gluon equation (19) are: 

R(k,q) = -± p U. (A.3) 

M T (k, q) = -2 ^ -^ ((k ■ q) 2 + k 2 f + 2f{k 2 + ?)) (A.4) 
k q p 

(qok-q-k f) Uk-p) 2 + k 2 f) 

Mi(k, q) = -2± L (A.5) 

(p k-p - kup 2 ) 2 ((k- q) 2 + k 2 c?) 

Mi(k, q) = -2A J V )_ (A . 6) 

(q^k 2 -(k-q) 2 )(q 2 p 2 -q p qpf 

M L (k, q) = -2± - — L . (A.7) 

k 2 q^p 1 k 2 p 2 
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The kernels for the heat-bath longitudinal part (20) have the form: 



P(k,q 

N T {k,q 

Ni(k,q 

N 2 (k,q 
N L (k,q 



= -2 



(q k 2 - k k ■ qj 
k 2 P 

(k k ■ q - q k 2 ^j ((q '■ p) 2 + cfp 2 ) 

(fk 2 — (k ■ q) 2 ^ (k 2 q 2 — k q kq) 2 

k 2 q 1 p 1 k 2 p 2 
q^k 2 -(k-q)A (k 2 p 2 - k oPo k P y 



k^cpfPktq 2 

(p-k)f (k p 2 +p k'' 



2k 2 q 2 p 2 k 2 q 2 p 2 
(p-q)k 2 (p q 2 + q p 2 ) + (k ■ q)f (k q 2 - q k 2 ^j 



(A.8) 
(A.9) 
(A.10) 
(A.ll) 

(A.12) 



B Numerical method 



B.l The algorithm 



The numerical method is based on approximating the unknown exact solution 
of the coupled system of integral equations by an iterative procedure. Fixed- 
point iteration works fine enough for simple cases, but in the present example 
suffers from extremely poor convergence. We therefore employed the Newton- 
Raphson method [37], which improves the convergence drastically, see e.g. 
ref. [7] (where this method has been employed for the first time for DSEs in 
the continuum) and ref. [16]. In the symbolic notation of eq. (21) we can write 
the DSEs as 



/(A 



l 



-U D (k 2 ) + U D (fi 2 ) = 0, 



(B.l) 



with Di = D(k{). The Newton-Raphson method then relies on linearly im- 
proved iteration steps 



£>r 1 = m - sDr\ sDr 1 = r\mfm 



(B.2) 



with 



J ( D i) 



dim 



(B.3) 
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We employ the following convergence criterion: 



< e for i e 1 . . . N . (B.4) 

Note that one should not use a convergence criterion for an averaged relative 
change in the variables as the dressing functions D(kf) are rather steep in the 
infrared, where only a few points ki are evaluated. Thus one can have situations 
where the majority of variables Di with ki in the ultraviolet are converged, 
whereas those A which correspond to small ki still change drastically. 

Furthermore, we observed that the algorithm converges much better when the 
Newton step is only performed partly when one is still far away from the exact 
solution. The weight c^Ncwton is calculated as 

ie{l,.../V}} . (B.5) 

If c^Newton > 1 it is set to one. This condition also ensures that the dressing 
functions are positive definite on all intermediate iteration steps. 

We determine the inverse of the matrix J(D^) approximately by employing 
a generalisation of the Gauss-Seidel method called 'successive overrelaxation' 
(SOR). This procedure is much faster than the exact inversion of the matrix. 
An exact determination of J{D 1 l)~ 1 is not needed anyway, since the inverse 
determines just the next respective iteration step, but is not used in calculating 
the right hand side of the integral equations. 

The SOR-algorithm necessarily converges, if the matrix J is strictly diagonally 
dominated, i.e. 

N 

\Ju\ >e + El J d ■ ( B - 6 ) 

The algorithm starts with an initial guess for the solution of the linear equation 
Jx = b with x, b e R N and J E R NxN . This guess is improved at each SOR- 
step 77 — > 77 + I according to: 

= (1 - *K + T [bi ~ E J^T 1 - E • (B.7) 

J%% \ J'=l j=i+l / 

The relaxation parameter uj interpolates linearly between the old vector x v 
and the Gauss-Seidel rule ou / ' Ju{- • •). The algorithm converges for values of 
u in the range < uj < 2. If the matrix J is strictly diagonally dominant, 
then convergence is best for 1 < u < 2 (i.e. 'overrelaxation'). However, other 
cases may give best results for a value between < uj < 1 ('underrelaxation'). 




^Newton = 0.99 Him 



5D 



r;+l 
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We tried values for u between 0.7 and 1.2 and obtained best results for w ~ 
0.8. This clearly indicates, that the Jacobian matrix is not strictly diagonally 
dominant but still well behaved enough for the SOR algorithm to converge. 
For the update of the i-th component of x, its relative change is 

/ i-l N 

V +1 = r [bi ~ E J^T 1 ~ E J^l I • (B.8) 

°i \ j= i j =i 

These quantities are used to calculate an estimated averaged error from the 
sum of the absolute values of the relative changes in the current approximation 
to the solution 

1 N , 

^ +1 = ^E|^ +1 | • (B.9) 
i=i 

The criterion for convergence is, that e v+1 < e, for a given small e. 

In our implementation of the Newton-Raphson method the SOR-algorithm 
has to be called at every step of the Newton-iteration. In general it delivers 
a converged result for J -1 . In the rare occasions where the SOR-algorithm 
did not converge, we determined J -1 with a Gauss- Jordan elimination and 
returned to the SOR-algorithm in the next Newton-step. 



B.2 Start guess for the dressing functions 



The Newton-Raphson iteration method described in the last section is a lo- 
cal method, which converges quadratically, if the start guess for the dressing 
functions is close enough to the final solutions. The exact meaning of 'close 
enough', however, depends very much on the problem at hand. If the dressing 
functions are infrared finite, as for example in the Dyson- Schwinger equations 
for the quark propagator [38,39], then convergence is achieved for almost any 
form of the start guess. However, if the dressing functions are infrared singu- 
lar, as in the present case, much more care is required and one should provide 
a start guess that is qualitatively similar to the actual solution 2 . In our case 
such a start guess can be provided by noting that the continuum version of the 
DSEs (18), (19) and (20) can be solved analytically in the zero temperature 
limit for both very large and very small momenta. 

In the infrared asymptotic limit, the ghost and gluon dressing functions, G{p 2 ) 
and Z(p 2 ), are proportional to simple power laws [6] 

Z( P 2 ) ~ ( P 2 ) 2K , G(p 2 ) ~ (p 2 )-, (B.10) 
2 As an alternative one can resort to a globally converging method, see ref. [16]. 
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Parameters 


C 


D s 


P 


V 


V 


G( P 2 ) 


0.35 


0.02 


-0.3 


0.5 


-0.2 


Z m (p 2 ),Z (p 2 ) 


550 


0.01 


1.2 


2.0 


-0.6 



Table B.l 

Parameters for the start guesses for the gluon and ghost dressing functions. 



93-V1201 
98 



0.595353 [5,9]. Correspondingly the running coupling 



with k 

(10) has an infrared fixed point, 

4tt r(3 - 2«)r(3 + «)r(i + «) 8.92 



a(0) 



QN r . 



r 2 (2 - «)r(2«) 



(B.ll) 



Both, the values of the exponent k and the fixed point a(0) have recently 
been shown to be independent of the gauge parameter in a class of gauges 
that interpolate between Landau and Coulomb gauge [32]. 

In the ultraviolet asymptotic limit one recovers the well known expressions 
from resummed perturbation theory, 



Z{p 2 ) = Z{p?) 
G{p 2 ) = G{n 2 ) 



llN c a(fi 2 ) 
12^ 

117V>(y) 

127T 



log ^ +1 



log 



' p 2 



+ 1 



(B.12) 
(B.13) 



Here Z(fi 2 ) and G(fi 2 ) denote the value of the dressing functions at some 
renormalisation point \j? and 7 and 5 are the respective anomalous dimensions. 
To one loop order one has 5 = —9/44 and 7 = — 1 — 25 = —13/22 for arbitrary 
number of colours N c and no quarks, Nf = 0. 



A simple ansatz that interpolates between both limits is given by 
F (P 2 ) = C ( , 7 ) f^log (p 2 /D s + l) + lV , 



(B.14) 



where C, D s , {3,v,rj are parameters to be chosen appropriately. Our choice of 
parameters for the ghost dressing function and the heat bath transverse and 
longitudinal dressing of the gluon are listed in table 3. As long as the qualita- 
tive properties of the start guess remains unchanged for the relevant momenta 
convergence is achieved for a whole range of parameter choices centred around 
the values given. 
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